ENM and SDM using ENMTools
Loading data
Just in case we don’t have our data saved from the last round:
library(dismo)## Loading required package: raster## Loading required package: splibrary(rgeos)## rgeos version: 0.5-3, (SVN revision 634)
## GEOS runtime version: 3.7.2-CAPI-1.11.2
## Linking to sp version: 1.4-1
## Polygon checking: TRUElibrary(rgbif)
library(leaflet)
ibm <- occ_search(scientificName = "Iberolacerta monticola", limit = 1500)
keep.cols <- c("species", "decimalLatitude", "decimalLongitude")
ibm <- ibm$data[,keep.cols]
ibm <- as.data.frame(unique(ibm))
ibm <- ibm[complete.cases(ibm),]
colnames(ibm) <- c("species", "lat", "lon")
head(ibm)## species lat lon
## 1 Iberolacerta monticola 43.45172 -6.699292
## 2 Iberolacerta monticola 43.43051 -7.188726
## 3 Iberolacerta monticola 43.46805 -6.703289
## 4 Iberolacerta monticola 43.51037 -6.504768
## 5 Iberolacerta monticola 43.52801 -6.556106
## 6 Iberolacerta monticola 43.55888 -6.728039all.worldclim <- raster::getData("worldclim", res = 10, var = "bio")
spain.worldclim <- crop(all.worldclim, extent(-10, 4, 35, 45))
m <- leaflet(ibm) %>%
addProviderTiles(provider = "Stamen.TerrainBackground") %>%
addRasterImage(spain.worldclim[[1]]) %>%
addCircleMarkers(lng = ~lon, lat = ~lat, radius = 5, stroke = FALSE,
fillOpacity = 1)
mOops! We’ve got one point down there on the west coast of Africa. We could explore why that is, or we could just eliminate it. For the sake of brevity, we’ll do the latter.
ibm <- ibm %>%
dplyr::filter(lat > 30)
m <- leaflet(ibm) %>%
addProviderTiles(provider = "Stamen.TerrainBackground") %>%
addRasterImage(spain.worldclim[[1]]) %>%
addCircleMarkers(lng = ~lon, lat = ~lat, radius = 5, stroke = FALSE,
fillOpacity = 1)
mNice!
Loading ENMTools
Now we’ll load ENMTools. If you’ve already installed it and run install.extras(), you can run it just by typing:
library(ENMTools)If you haven’t done either of those things, see the “Installing ENMTools” tutorial for details.
Building an enmtools.species object
ENMTools automates a lot of the fiddly parts of building SDMs, but in order to do that it requires us to get our data into a specific format.
monticola <- enmtools.species()
monticola##
##
## Range raster not defined.
##
## Presence points not defined.
##
## Background points not defined.
##
## Species name not defined.We now have an empty enmtools.species object. It’s basically just a list where we can store our data. This makes sure that our data is in a predictable format for ENMTools’ modeling functions. Let’s go ahead and add data to it now.
monticola$presence.points <- ibm[,c("lon", "lat")]
monticola$species.name <- "Iberolacerta monticola"
monticola##
##
## Range raster not defined.
##
## Presence points (first ten only):
##
## lon lat
## ---------- ---------
## -6.699292 43.45172
## -7.188726 43.43051
## -6.703289 43.46805
## -6.504768 43.51037
## -6.556106 43.52801
## -6.728039 43.55888
## -6.823801 43.39811
## -6.732560 43.54083
## -6.650406 43.40437
## -6.541960 43.45676
##
##
## Background points not defined.
##
## Species name: Iberolacerta monticolaWe can define a range for our species using any raster, or we can build one using the background.raster.buffer function of ENMTools
monticola$range <- background.raster.buffer(monticola$presence.points, 50000, mask = spain.worldclim)
plot(monticola)
This range will be used for sampling background points in ENM construction, but we could also provide points manually if we wanted to using monticola$background.points <- something.
It’s always a good idea to run your species through the “check.species” function to make sure it’s formatted correctly
monticola <- check.species(monticola)Not only does this tell us if we’ve messed something up, it reformats column names etc. so that they’re in the appropriate format for the remaining ENMTools features. The fact that everything is stored in a reliable place in these enmtools.species objects allows us to streamline a lot of operations. For instance, instead of manually building a leaflet plot we can now just use the enmtools.species object to build one automatically.
interactive.plot.enmtools.species(monticola)Building ENMs with ENMTools
Okay, we’ve got a species. Let’s build a model! We’ll start with the two most popular climate envelope models.
monticola.bc <- enmtools.bc(species = monticola, env = spain.worldclim, test.prop = 0.3, nback = 400)##
##
## No background points provided, drawing background from range raster.## Warning in randomPoints(species$range, nback, species$presence.points):
## generated random points = 0.9975 times requested number##
##
## Drawing background from species background points.##
##
## Drawing background from species background points.monticola.bc##
##
## Data table (top ten lines):
##
## Longitude Latitude
## --- ---------- ---------
## 1 -6.699292 43.45172
## 2 -7.188726 43.43051
## 3 -6.703289 43.46805
## 6 -6.728039 43.55888
## 10 -6.541960 43.45676
## 11 -6.910225 43.29794
## 12 -6.659564 43.48790
## 13 -6.817537 43.28005
## 14 -6.831939 43.38955
## 15 -6.793121 43.45933
##
##
## Model: class : Bioclim
##
## variables: bio1 bio2 bio3 bio4 bio5 bio6 bio7 bio8 bio9 bio10 bio11 bio12 bio13 bio14 bio15 bio16 bio17 bio18 bio19
##
##
## presence points: 246
## bio1 bio2 bio3 bio4 bio5 bio6 bio7 bio8 bio9 bio10 bio11 bio12 bio13 bio14
## 1 121 88 41 4237 239 28 211 97 175 178 70 920 121 36
## 2 128 84 41 4080 239 38 201 83 180 182 78 959 124 33
## 3 121 88 41 4237 239 28 211 97 175 178 70 920 121 36
## 4 137 83 42 3957 244 48 196 115 187 190 89 852 114 33
## 5 117 90 42 4292 237 23 214 93 172 174 65 926 123 38
## 6 117 91 41 4472 240 20 220 67 175 176 63 942 122 34
## 7 117 90 42 4292 237 23 214 93 172 174 65 926 123 38
## 8 105 93 41 4614 233 7 226 54 165 166 49 983 127 38
## 9 121 88 41 4237 239 28 211 97 175 178 70 920 121 36
## 10 121 88 41 4237 239 28 211 97 175 178 70 920 121 36
## bio15 bio16 bio17 bio18 bio19
## 1 31 320 135 148 279
## 2 34 345 126 142 310
## 3 31 320 135 148 279
## 4 31 298 124 139 253
## 5 30 320 140 153 274
## 6 33 334 129 142 300
## 7 30 320 140 153 274
## 8 32 345 141 152 309
## 9 31 320 135 148 279
## 10 31 320 135 148 279
## (... ... ...)
##
##
##
## Model fit (training data): class : ModelEvaluation
## n presences : 246
## n absences : 399
## AUC : 0.8480958
## cor : 0.5188009
## max TPR+TNR at : 0.04868049
##
##
## Environment space model fit (training data): class : ModelEvaluation
## n presences : 246
## n absences : 10000
## AUC : 0.9995441
## cor : 0.7325
## max TPR+TNR at : 0.00396504
##
##
## Proportion of data wittheld for model fitting: 0.3
##
## Model fit (test data): class : ModelEvaluation
## n presences : 106
## n absences : 399
## AUC : 0.8844162
## cor : 0.5681143
## max TPR+TNR at : 0.07713577
##
##
## Environment space model fit (test data): class : ModelEvaluation
## n presences : 106
## n absences : 10000
## AUC : 0.9998014
## cor : 0.7681282
## max TPR+TNR at : 0.00803008
##
##
## Suitability:
## class : RasterLayer
## dimensions : 60, 84, 5040 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -10, 4, 35, 45 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : layer
## values : 0, 0.601626 (min, max)
##
##
##
## Notes:
## NULL
monticola.bc$response.plots[3]## $bio3
monticola.bc$response.plots[8]## $bio8
visualize.enm(monticola.bc, spain.worldclim, layers = c("bio3", "bio8"))## $background.plot
## [1] NA
##
## $suit.plot
monticola.dm <- enmtools.dm(species = monticola, env = spain.worldclim, test.prop = 0.3, nback = 400)##
##
## No background points provided, drawing background from range raster.## Warning in randomPoints(species$range, nback, species$presence.points):
## generated random points = 0.9975 times requested number##
##
## Drawing background from species background points.##
##
## Drawing background from species background points.monticola.dm##
##
## Data table (top ten lines):
##
## Longitude Latitude
## --- ---------- ---------
## 1 -6.699292 43.45172
## 2 -7.188726 43.43051
## 3 -6.703289 43.46805
## 4 -6.504768 43.51037
## 5 -6.556106 43.52801
## 6 -6.728039 43.55888
## 7 -6.823801 43.39811
## 8 -6.732560 43.54083
## 9 -6.650406 43.40437
## 11 -6.910225 43.29794
##
##
## Model: class : Domain
##
## variables: bio1 bio2 bio3 bio4 bio5 bio6 bio7 bio8 bio9 bio10 bio11 bio12 bio13 bio14 bio15 bio16 bio17 bio18 bio19
##
##
## presence points: 246
## bio1 bio2 bio3 bio4 bio5 bio6 bio7 bio8 bio9 bio10 bio11 bio12 bio13 bio14
## 1 121 88 41 4237 239 28 211 97 175 178 70 920 121 36
## 2 128 84 41 4080 239 38 201 83 180 182 78 959 124 33
## 3 121 88 41 4237 239 28 211 97 175 178 70 920 121 36
## 4 134 85 42 3987 243 45 198 113 185 187 86 860 116 35
## 5 134 85 42 3987 243 45 198 113 185 187 86 860 116 35
## 6 137 83 42 3957 244 48 196 115 187 190 89 852 114 33
## 7 121 88 41 4237 239 28 211 97 175 178 70 920 121 36
## 8 137 83 42 3957 244 48 196 115 187 190 89 852 114 33
## 9 117 90 42 4292 237 23 214 93 172 174 65 926 123 38
## 10 117 91 41 4472 240 20 220 67 175 176 63 942 122 34
## bio15 bio16 bio17 bio18 bio19
## 1 31 320 135 148 279
## 2 34 345 126 142 310
## 3 31 320 135 148 279
## 4 30 299 130 145 248
## 5 30 299 130 145 248
## 6 31 298 124 139 253
## 7 31 320 135 148 279
## 8 31 298 124 139 253
## 9 30 320 140 153 274
## 10 33 334 129 142 300
## (... ... ...)
##
##
##
## Model fit (training data): class : ModelEvaluation
## n presences : 246
## n absences : 399
## AUC : 0.7956782
## cor : 0.4647681
## max TPR+TNR at : 0.5332631
##
##
## Environment space model fit (training data): class : ModelEvaluation
## n presences : 246
## n absences : 10000
## AUC : 0.9973431
## cor : 0.4724289
## max TPR+TNR at : 0.4364953
##
##
## Proportion of data wittheld for model fitting: 0.3
##
## Model fit (test data): class : ModelEvaluation
## n presences : 106
## n absences : 399
## AUC : 0.7857025
## cor : 0.3637048
## max TPR+TNR at : 0.6109946
##
##
## Environment space model fit (test data): class : ModelEvaluation
## n presences : 106
## n absences : 10000
## AUC : 0.9921604
## cor : 0.33623
## max TPR+TNR at : 0.4485913
##
##
## Suitability:
## class : RasterLayer
## dimensions : 60, 84, 5040 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -10, 4, 35, 45 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : layer
## values : 0, 0.771146 (min, max)
##
##
##
## Notes:
## NULL
monticola.dm$response.plots[3]## $bio3
monticola.dm$response.plots[8]## $bio8
visualize.enm(monticola.dm, spain.worldclim, layers = c("bio3", "bio8"))## $background.plot
## [1] NA
##
## $suit.plot
monticola.glm <- enmtools.glm(species = monticola, env = spain.worldclim)##
##
## No background points provided, drawing background from range raster.## Adding environmental data to species Iberolacerta monticola## Processing presence points...## Processing background points...## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
## prediction from a rank-deficient fit may be misleading
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
## prediction from a rank-deficient fit may be misleading
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
## prediction from a rank-deficient fit may be misleading
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
## prediction from a rank-deficient fit may be misleading
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
## prediction from a rank-deficient fit may be misleading
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
## prediction from a rank-deficient fit may be misleading##
##
## Drawing background from species background points.## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
## prediction from a rank-deficient fit may be misleading
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
## prediction from a rank-deficient fit may be misleading
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
## prediction from a rank-deficient fit may be misleading
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
## prediction from a rank-deficient fit may be misleading
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
## prediction from a rank-deficient fit may be misleading
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
## prediction from a rank-deficient fit may be misleading
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
## prediction from a rank-deficient fit may be misleading
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
## prediction from a rank-deficient fit may be misleading
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
## prediction from a rank-deficient fit may be misleading
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
## prediction from a rank-deficient fit may be misleading
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
## prediction from a rank-deficient fit may be misleading
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
## prediction from a rank-deficient fit may be misleading
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
## prediction from a rank-deficient fit may be misleading
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
## prediction from a rank-deficient fit may be misleading
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
## prediction from a rank-deficient fit may be misleading
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
## prediction from a rank-deficient fit may be misleading
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
## prediction from a rank-deficient fit may be misleading
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
## prediction from a rank-deficient fit may be misleading
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
## prediction from a rank-deficient fit may be misleading
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
## prediction from a rank-deficient fit may be misleading
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
## prediction from a rank-deficient fit may be misleadingwarnings()We’ve got a warning here, which is due to the same reason as before - too many predictors and not enough data.
Let’s look at the model first and then we’ll work on that.
monticola.glm##
##
## Formula: presence ~ bio1 + bio2 + bio3 + bio4 + bio5 + bio6 + bio7 + bio8 +
## bio9 + bio10 + bio11 + bio12 + bio13 + bio14 + bio15 + bio16 +
## bio17 + bio18 + bio19
## <environment: 0x7fcc5746dac8>
##
##
## Data table (top ten lines):
##
## Longitude Latitude bio1 bio2 bio3 bio4 bio5 bio6 bio7 bio8 bio9 bio10 bio11 bio12 bio13 bio14 bio15 bio16 bio17 bio18 bio19 presence
## ---------- --------- ----- ----- ----- ----- ----- ----- ----- ----- ----- ------ ------ ------ ------ ------ ------ ------ ------ ------ ------ ---------
## -6.699292 43.45172 121 88 41 4237 239 28 211 97 175 178 70 920 121 36 31 320 135 148 279 1
## -7.188726 43.43051 128 84 41 4080 239 38 201 83 180 182 78 959 124 33 34 345 126 142 310 1
## -6.703289 43.46805 121 88 41 4237 239 28 211 97 175 178 70 920 121 36 31 320 135 148 279 1
## -6.504768 43.51037 134 85 42 3987 243 45 198 113 185 187 86 860 116 35 30 299 130 145 248 1
## -6.556106 43.52801 134 85 42 3987 243 45 198 113 185 187 86 860 116 35 30 299 130 145 248 1
## -6.728039 43.55888 137 83 42 3957 244 48 196 115 187 190 89 852 114 33 31 298 124 139 253 1
## -6.823801 43.39811 121 88 41 4237 239 28 211 97 175 178 70 920 121 36 31 320 135 148 279 1
## -6.732560 43.54083 137 83 42 3957 244 48 196 115 187 190 89 852 114 33 31 298 124 139 253 1
## -6.650406 43.40437 117 90 42 4292 237 23 214 93 172 174 65 926 123 38 30 320 140 153 274 1
## -6.541960 43.45676 117 90 42 4292 237 23 214 93 172 174 65 926 123 38 30 320 140 153 274 1
##
##
## Model:
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1,
## 2)], weights = weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7585 -0.7792 -0.2342 0.7474 3.2781
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.486842 14.349989 -0.452 0.651236
## bio1 -0.655210 0.186813 -3.507 0.000453 ***
## bio2 0.272431 0.163512 1.666 0.095689 .
## bio3 0.326360 0.291146 1.121 0.262310
## bio4 0.012902 0.004887 2.640 0.008291 **
## bio5 -0.280712 0.100930 -2.781 0.005415 **
## bio6 0.092671 0.109073 0.850 0.395537
## bio7 NA NA NA NA
## bio8 -0.022938 0.008745 -2.623 0.008714 **
## bio9 0.035934 0.067565 0.532 0.594831
## bio10 0.038312 0.223170 0.172 0.863697
## bio11 0.740189 0.221573 3.341 0.000836 ***
## bio12 -0.010986 0.015142 -0.726 0.468112
## bio13 -0.029120 0.042102 -0.692 0.489164
## bio14 -0.221261 0.151081 -1.465 0.143053
## bio15 -0.167661 0.079782 -2.101 0.035598 *
## bio16 0.085904 0.031099 2.762 0.005740 **
## bio17 0.016308 0.067458 0.242 0.808976
## bio18 0.019858 0.037005 0.537 0.591527
## bio19 -0.035383 0.028793 -1.229 0.219105
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 975.95 on 870 degrees of freedom
## Residual deviance: 711.76 on 852 degrees of freedom
## AIC: 924.16
##
## Number of Fisher Scoring iterations: 8
##
##
##
## Model fit (training data): class : ModelEvaluation
## n presences : 352
## n absences : 519
## AUC : 0.8258561
## cor : 0.5013653
## max TPR+TNR at : -0.05890047
##
##
## Environment space model fit (training data): class : ModelEvaluation
## n presences : 352
## n absences : 10000
## AUC : 0.5502134
## cor : 0.07998986
## max TPR+TNR at : 0.1526213
##
##
## Proportion of data wittheld for model fitting: 0
##
## Model fit (test data): [1] NA
##
##
## Environment space model fit (test data): [1] NA
##
##
## Suitability:
## class : RasterLayer
## dimensions : 60, 84, 5040 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -10, 4, 35, 45 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : layer
## values : 4.14525e-06, 0.9960972 (min, max)
##
##
##
## Notes:
names(monticola.glm)## [1] "species.name" "formula"
## [3] "analysis.df" "test.data"
## [5] "test.prop" "model"
## [7] "training.evaluation" "test.evaluation"
## [9] "env.training.evaluation" "env.test.evaluation"
## [11] "rts.test" "suitability"
## [13] "call" "notes"
## [15] "response.plots"Look at all of that stuff!
monticola.glm$response.plots## $bio1
##
## $bio2
##
## $bio3
##
## $bio4
##
## $bio5
##
## $bio6
##
## $bio7
##
## $bio8
##
## $bio9
##
## $bio10
##
## $bio11
##
## $bio12
##
## $bio13
##
## $bio14
##
## $bio15
##
## $bio16
##
## $bio17
##
## $bio18
##
## $bio19
We have interactive plots for models too!
interactive.plot.enmtools.model(monticola.glm)Let’s try fitting a simpler model
First let’s look at the correlations between variables
raster.cor.matrix(spain.worldclim)## bio1 bio2 bio3 bio4 bio5 bio6
## bio1 1.000000000 -0.05192243 -0.02536989 -0.06482844 0.7309560 0.89386921
## bio2 -0.051922434 1.00000000 0.25804128 0.68346008 0.5595536 -0.41646754
## bio3 -0.025369890 0.25804128 1.00000000 -0.50487533 -0.1563297 0.09040494
## bio4 -0.064828440 0.68346008 -0.50487533 1.00000000 0.5735377 -0.47341350
## bio5 0.730955951 0.55955362 -0.15632969 0.57353766 1.0000000 0.37935531
## bio6 0.893869213 -0.41646754 0.09040494 -0.47341350 0.3793553 1.00000000
## bio7 -0.031962836 0.88012003 -0.22456546 0.94036051 0.6494830 -0.45715397
## bio8 0.404382488 -0.11713369 -0.09355248 -0.03873301 0.2024989 0.33552051
## bio9 0.683287223 0.20191286 0.01963002 0.11583860 0.6773847 0.54402099
## bio10 0.922563848 0.20792847 -0.22821990 0.32228650 0.9151040 0.66801013
## bio11 0.946033848 -0.27710014 0.11951392 -0.37865009 0.4946439 0.98216169
## bio12 -0.392848782 -0.45597177 0.25444757 -0.57123600 -0.6674771 -0.06869193
## bio13 -0.120028527 -0.51403698 0.25059275 -0.64521027 -0.4924741 0.20166588
## bio14 -0.758765734 -0.25110222 0.10775308 -0.24002532 -0.8123826 -0.56244009
## bio15 0.773143662 -0.07304360 -0.01876198 -0.12025929 0.5822825 0.74503621
## bio16 -0.087315073 -0.49267982 0.26595406 -0.63962927 -0.4507268 0.23058503
## bio17 -0.739725188 -0.26496834 0.12139640 -0.26115476 -0.8098741 -0.53463352
## bio18 -0.713437527 -0.33085626 0.08080540 -0.28904180 -0.8266616 -0.49637460
## bio19 -0.002591276 -0.40575986 0.29614464 -0.58921379 -0.3270008 0.28232294
## bio7 bio8 bio9 bio10 bio11 bio12
## bio1 -0.03196284 0.40438249 0.68328722 0.9225638 0.94603385 -0.39284878
## bio2 0.88012003 -0.11713369 0.20191286 0.2079285 -0.27710014 -0.45597177
## bio3 -0.22456546 -0.09355248 0.01963002 -0.2282199 0.11951392 0.25444757
## bio4 0.94036051 -0.03873301 0.11583860 0.3222865 -0.37865009 -0.57123600
## bio5 0.64948302 0.20249891 0.67738465 0.9151040 0.49464394 -0.66747713
## bio6 -0.45715397 0.33552051 0.54402099 0.6680101 0.98216169 -0.06869193
## bio7 1.00000000 -0.08108260 0.20404940 0.3306597 -0.33167441 -0.58515377
## bio8 -0.08108260 1.00000000 -0.04153511 0.3626530 0.38318128 -0.42603645
## bio9 0.20404940 -0.04153511 1.00000000 0.7008436 0.60235162 -0.31651363
## bio10 0.33065966 0.36265299 0.70084359 1.0000000 0.75367532 -0.59770879
## bio11 -0.33167441 0.38318128 0.60235162 0.7536753 1.00000000 -0.18778961
## bio12 -0.58515377 -0.42603645 -0.31651363 -0.5977088 -0.18778961 1.00000000
## bio13 -0.63911541 -0.33110726 -0.12411293 -0.3619468 0.09551341 0.92211902
## bio14 -0.31867777 -0.22120096 -0.66922270 -0.8260330 -0.64267656 0.62160840
## bio15 -0.05256181 0.06324795 0.60151828 0.7025760 0.77296096 -0.20402064
## bio16 -0.62275227 -0.37105106 -0.08385793 -0.3281109 0.12411688 0.91197261
## bio17 -0.33911807 -0.22760761 -0.63176811 -0.8165397 -0.61836117 0.66182192
## bio18 -0.38669626 -0.13863713 -0.68989269 -0.8016505 -0.58312240 0.64972769
## bio19 -0.54634045 -0.44956955 0.05633320 -0.2288494 0.18398498 0.85972346
## bio13 bio14 bio15 bio16 bio17 bio18
## bio1 -0.12002853 -0.7587657 0.77314366 -0.08731507 -0.7397252 -0.7134375
## bio2 -0.51403698 -0.2511022 -0.07304360 -0.49267982 -0.2649683 -0.3308563
## bio3 0.25059275 0.1077531 -0.01876198 0.26595406 0.1213964 0.0808054
## bio4 -0.64521027 -0.2400253 -0.12025929 -0.63962927 -0.2611548 -0.2890418
## bio5 -0.49247406 -0.8123826 0.58228255 -0.45072685 -0.8098741 -0.8266616
## bio6 0.20166588 -0.5624401 0.74503621 0.23058503 -0.5346335 -0.4963746
## bio7 -0.63911541 -0.3186778 -0.05256181 -0.62275227 -0.3391181 -0.3866963
## bio8 -0.33110726 -0.2212010 0.06324795 -0.37105106 -0.2276076 -0.1386371
## bio9 -0.12411293 -0.6692227 0.60151828 -0.08385793 -0.6317681 -0.6898927
## bio10 -0.36194679 -0.8260330 0.70257598 -0.32811094 -0.8165397 -0.8016505
## bio11 0.09551341 -0.6426766 0.77296096 0.12411688 -0.6183612 -0.5831224
## bio12 0.92211902 0.6216084 -0.20402064 0.91197261 0.6618219 0.6497277
## bio13 1.00000000 0.3184485 0.13929966 0.99092613 0.3638003 0.3707213
## bio14 0.31844853 1.0000000 -0.83281349 0.27096316 0.9913040 0.9791848
## bio15 0.13929966 -0.8328135 1.00000000 0.19326285 -0.8288811 -0.8093866
## bio16 0.99092613 0.2709632 0.19326285 1.00000000 0.3156585 0.3166368
## bio17 0.36380031 0.9913040 -0.82888105 0.31565850 1.0000000 0.9832066
## bio18 0.37072125 0.9791848 -0.80938657 0.31663682 0.9832066 1.0000000
## bio19 0.94856030 0.1661662 0.28320125 0.97382680 0.2114466 0.1873645
## bio19
## bio1 -0.002591276
## bio2 -0.405759856
## bio3 0.296144638
## bio4 -0.589213793
## bio5 -0.327000835
## bio6 0.282322938
## bio7 -0.546340454
## bio8 -0.449569548
## bio9 0.056333201
## bio10 -0.228849357
## bio11 0.183984983
## bio12 0.859723460
## bio13 0.948560299
## bio14 0.166166182
## bio15 0.283201249
## bio16 0.973826798
## bio17 0.211446620
## bio18 0.187364463
## bio19 1.000000000That’s handy, but maybe not as intuitive as plotting it
raster.cor.plot(spain.worldclim)## $cor.mds.plot
##
## $cor.heatmap
This contains both a heatmap and an mds plot, let’s look at both!
Okay, now a simpler GLM.
monticola.glm <- enmtools.glm(monticola, spain.worldclim, f = pres ~ poly(bio3, 2) + poly(bio8, 2) + poly(bio11, 2), nback = 400)##
##
## No background points provided, drawing background from range raster.## Warning in randomPoints(species$range, nback, species$presence.points):
## generated random points = 0.9975 times requested number## Adding environmental data to species Iberolacerta monticola## Processing presence points...## Processing background points...##
##
## Drawing background from species background points.interactive.plot.enmtools.model(monticola.glm)monticola.glm$response.plots## $bio3
##
## $bio8
##
## $bio11
Notice that we have flat lines for most of the responses - that’s because we only built a model using variables 1 and 8, but passed it all 19.
Let’s look at what our model looks like in environment space!
visualize.enm(model = monticola.glm, env = spain.worldclim, layers = c("bio8", "bio11"))## $background.plot## Warning: Removed 396 rows containing missing values (geom_raster).
##
## $suit.plot
Now let’s try one with some test data for evaluation
monticola.glm <- enmtools.glm(monticola, spain.worldclim, f = pres ~ poly(bio3, 2) + poly(bio8, 2) + poly(bio11, 2), nback = 400, test.prop = 0.3)##
##
## No background points provided, drawing background from range raster.## Warning in randomPoints(species$range, nback, species$presence.points):
## generated random points = 0.9975 times requested number## Adding environmental data to species Iberolacerta monticola## Processing presence points...## Processing background points...##
##
## Drawing background from species background points.
##
##
##
##
## Drawing background from species background points.monticola.glm##
##
## Formula: presence ~ poly(bio3, 2) + poly(bio8, 2) + poly(bio11, 2)
## <environment: 0x7fcc6946b070>
##
##
## Data table (top ten lines):
##
## Longitude Latitude bio1 bio2 bio3 bio4 bio5 bio6 bio7 bio8 bio9 bio10 bio11 bio12 bio13 bio14 bio15 bio16 bio17 bio18 bio19 presence
## --- ---------- --------- ----- ----- ----- ----- ----- ----- ----- ----- ----- ------ ------ ------ ------ ------ ------ ------ ------ ------ ------ ---------
## 4 -6.504768 43.51037 134 85 42 3987 243 45 198 113 185 187 86 860 116 35 30 299 130 145 248 1
## 5 -6.556106 43.52801 134 85 42 3987 243 45 198 113 185 187 86 860 116 35 30 299 130 145 248 1
## 6 -6.728039 43.55888 137 83 42 3957 244 48 196 115 187 190 89 852 114 33 31 298 124 139 253 1
## 8 -6.732560 43.54083 137 83 42 3957 244 48 196 115 187 190 89 852 114 33 31 298 124 139 253 1
## 10 -6.541960 43.45676 117 90 42 4292 237 23 214 93 172 174 65 926 123 38 30 320 140 153 274 1
## 13 -6.817537 43.28005 105 93 41 4614 233 7 226 54 165 166 49 983 127 38 32 345 141 152 309 1
## 14 -6.831939 43.38955 121 88 41 4237 239 28 211 97 175 178 70 920 121 36 31 320 135 148 279 1
## 17 -6.619466 43.42050 117 90 42 4292 237 23 214 93 172 174 65 926 123 38 30 320 140 153 274 1
## 18 -6.724369 43.52515 137 83 42 3957 244 48 196 115 187 190 89 852 114 33 31 298 124 139 253 1
## 19 -7.657592 40.33878 101 76 33 5359 236 7 229 39 173 173 39 1514 234 15 57 621 92 92 621 1
##
##
## Model:
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1,
## 2)], weights = weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3298 -0.6546 -0.3187 0.6020 2.5617
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.3104 0.1231 -2.522 0.01168 *
## poly(bio3, 2)1 37.1228 3.8332 9.684 < 2e-16 ***
## poly(bio3, 2)2 22.1697 3.7077 5.979 2.24e-09 ***
## poly(bio8, 2)1 -21.8758 4.0638 -5.383 7.32e-08 ***
## poly(bio8, 2)2 9.7790 3.3682 2.903 0.00369 **
## poly(bio11, 2)1 -6.8819 4.6121 -1.492 0.13566
## poly(bio11, 2)2 -9.8798 3.7171 -2.658 0.00786 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 682.06 on 644 degrees of freedom
## Residual deviance: 462.41 on 638 degrees of freedom
## AIC: 621.59
##
## Number of Fisher Scoring iterations: 5
##
##
##
## Model fit (training data): class : ModelEvaluation
## n presences : 246
## n absences : 399
## AUC : 0.8682529
## cor : 0.5992082
## max TPR+TNR at : 0.5106094
##
##
## Environment space model fit (training data): class : ModelEvaluation
## n presences : 246
## n absences : 10000
## AUC : 0.6365004
## cor : 0.09199756
## max TPR+TNR at : 0.4073416
##
##
## Proportion of data wittheld for model fitting: 0.3
##
## Model fit (test data): class : ModelEvaluation
## n presences : 106
## n absences : 399
## AUC : 0.8457937
## cor : 0.4924478
## max TPR+TNR at : 0.4290774
##
##
## Environment space model fit (test data): class : ModelEvaluation
## n presences : 106
## n absences : 10000
## AUC : 0.7015113
## cor : 0.07935285
## max TPR+TNR at : 0.6055772
##
##
## Suitability:
## class : RasterLayer
## dimensions : 60, 84, 5040 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -10, 4, 35, 45 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : layer
## values : 0.003412865, 1 (min, max)
##
##
##
## Notes:
interactive.plot.enmtools.model(monticola.glm)Hey wait! What are these environment space model fit metrics? We’ll get to that in a little while.
ENMTools also allows us to build models using other methods.
monticola.gam <- enmtools.gam(species = monticola, env = spain.worldclim, f = pres ~ bio3 + bio8 + bio11, test.prop = 0.3, nback = 400)##
##
## No background points provided, drawing background from range raster.## Warning in randomPoints(species$range, nback, species$presence.points):
## generated random points = 0.9975 times requested number## Adding environmental data to species Iberolacerta monticola## Processing presence points...## Processing background points...##
##
## Drawing background from species background points.
##
##
##
##
## Drawing background from species background points.monticola.gam##
##
## Formula: presence ~ bio3 + bio8 + bio11
## <environment: 0x7fcc50d850b8>
##
##
## Data table (top ten lines):
##
## Longitude Latitude bio1 bio2 bio3 bio4 bio5 bio6 bio7 bio8 bio9 bio10 bio11 bio12 bio13 bio14 bio15 bio16 bio17 bio18 bio19 presence
## --- ---------- --------- ----- ----- ----- ----- ----- ----- ----- ----- ----- ------ ------ ------ ------ ------ ------ ------ ------ ------ ------ ---------
## 3 -6.703289 43.46805 121 88 41 4237 239 28 211 97 175 178 70 920 121 36 31 320 135 148 279 1
## 4 -6.504768 43.51037 134 85 42 3987 243 45 198 113 185 187 86 860 116 35 30 299 130 145 248 1
## 5 -6.556106 43.52801 134 85 42 3987 243 45 198 113 185 187 86 860 116 35 30 299 130 145 248 1
## 6 -6.728039 43.55888 137 83 42 3957 244 48 196 115 187 190 89 852 114 33 31 298 124 139 253 1
## 7 -6.823801 43.39811 121 88 41 4237 239 28 211 97 175 178 70 920 121 36 31 320 135 148 279 1
## 8 -6.732560 43.54083 137 83 42 3957 244 48 196 115 187 190 89 852 114 33 31 298 124 139 253 1
## 9 -6.650406 43.40437 117 90 42 4292 237 23 214 93 172 174 65 926 123 38 30 320 140 153 274 1
## 10 -6.541960 43.45676 117 90 42 4292 237 23 214 93 172 174 65 926 123 38 30 320 140 153 274 1
## 11 -6.910225 43.29794 117 91 41 4472 240 20 220 67 175 176 63 942 122 34 33 334 129 142 300 1
## 14 -6.831939 43.38955 121 88 41 4237 239 28 211 97 175 178 70 920 121 36 31 320 135 148 279 1
##
##
## Model:
## Family: binomial
## Link function: logit
##
## Formula:
## presence ~ bio3 + bio8 + bio11
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -22.506470 2.611045 -8.620 < 2e-16 ***
## bio3 0.637652 0.069526 9.171 < 2e-16 ***
## bio8 -0.025985 0.005292 -4.910 9.12e-07 ***
## bio11 -0.008863 0.005484 -1.616 0.106
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## R-sq.(adj) = 0.305 Deviance explained = 23.6%
## -REML = 272.41 Scale est. = 1 n = 645
##
##
## Model fit (training data): class : ModelEvaluation
## n presences : 246
## n absences : 399
## AUC : 0.8320955
## cor : 0.5149226
## max TPR+TNR at : 0.4368832
##
##
## Environment space model fit (training data): class : ModelEvaluation
## n presences : 246
## n absences : 10000
## AUC : 0.6478854
## cor : 0.08988055
## max TPR+TNR at : 0.5642361
##
##
## Proportion of data wittheld for model fitting: 0.3
##
## Model fit (test data): class : ModelEvaluation
## n presences : 106
## n absences : 399
## AUC : 0.8458883
## cor : 0.4483307
## max TPR+TNR at : 0.3387611
##
##
## Environment space model fit (test data): class : ModelEvaluation
## n presences : 106
## n absences : 10000
## AUC : 0.7260755
## cor : 0.08606299
## max TPR+TNR at : 0.5260684
##
##
## Suitability:
## class : RasterLayer
## dimensions : 60, 84, 5040 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -10, 4, 35, 45 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : layer
## values : 6.582324e-05, 0.9932084 (min, max)
##
##
##
## Notes:
## NULL
monticola.gam$response.plots## $bio3
##
## $bio8
##
## $bio11
visualize.enm(monticola.gam, spain.worldclim, layers = c("bio11", "bio8"))## $background.plot## Warning: Removed 396 rows containing missing values (geom_raster).
##
## $suit.plot
If you have maxent installed and setup with dismo, you can even build maxent models.
monticola.mx <- enmtools.maxent(species = monticola, env = spain.worldclim, test.prop = 0.3, nback = 400)##
##
## No background points provided, drawing background from range raster.## Warning in randomPoints(species$range, nback, species$presence.points):
## generated random points = 0.9975 times requested number##
##
## Drawing background from species background points.##
##
## Drawing background from species background points.monticola.mx##
##
## Data table (top ten lines):
##
## Longitude Latitude presence
## --- ---------- --------- ---------
## 1 -6.699292 43.45172 1
## 2 -7.188726 43.43051 1
## 3 -6.703289 43.46805 1
## 10 -6.541960 43.45676 1
## 12 -6.659564 43.48790 1
## 13 -6.817537 43.28005 1
## 15 -6.793121 43.45933 1
## 16 -6.499937 43.47067 1
## 18 -6.724369 43.52515 1
## 21 -6.933462 43.48396 1
##
##
## Model: Length Class Mode
## 1 MaxEnt S4
##
##
## Model fit (training data): class : ModelEvaluation
## n presences : 246
## n absences : 399
## AUC : 0.9474703
## cor : 0.7935504
## max TPR+TNR at : 0.419713
##
##
## Environment space model fit (training data): class : ModelEvaluation
## n presences : 246
## n absences : 10000
## AUC : 0.6758703
## cor : 0.1143143
## max TPR+TNR at : 0.409728
##
##
## Proportion of data wittheld for model fitting: 0.3
##
## Model fit (test data): class : ModelEvaluation
## n presences : 106
## n absences : 399
## AUC : 0.9257341
## cor : 0.7022422
## max TPR+TNR at : 0.5160748
##
##
## Environment space model fit (test data): class : ModelEvaluation
## n presences : 106
## n absences : 10000
## AUC : 0.6538226
## cor : 0.06436646
## max TPR+TNR at : 0.5804668
##
##
## Suitability:
## class : RasterLayer
## dimensions : 60, 84, 5040 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -10, 4, 35, 45 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : layer
## values : 0.0008448937, 0.9982985 (min, max)
##
##
##
## Notes:
## NULLmonticola.mx$response.plots## $bio1
##
## $bio2
##
## $bio3
##
## $bio4
##
## $bio5
##
## $bio6
##
## $bio7
##
## $bio8
##
## $bio9
##
## $bio10
##
## $bio11
##
## $bio12
##
## $bio13
##
## $bio14
##
## $bio15
##
## $bio16
##
## $bio17
##
## $bio18
##
## $bio19
visualize.enm(monticola.mx, spain.worldclim, layers = c("bio3", "bio8"))## $background.plot## Warning: Removed 396 rows containing missing values (geom_raster).
##
## $suit.plot
Note that what ENMTools is doing for ENM construction is just acting as a front end for dismo and ppmlasso: the interface automates a bunch of stuff, but the underlying math is just dismo and ppmlasso. Stored in the enmtools.model object (monticola.mx), there is a model we can project using dismo’s predict function. We can even use it to bring up the maxent model page in a web browser.
You can even build poisson point process models using ppmlasso. Note: ppmlasso can take a LONG time to fit models with a bunch of predictors, so it’s best to pass a formula here instead of using all 19 variables.
monticola.ppm <- enmtools.ppmlasso(species = monticola, env = spain.worldclim, f = pres ~ bio3 + bio8 + bio11, test.prop = 0.3, nback = 400, gamma = 0)##
##
## No background points provided, drawing background from range raster.## Warning in randomPoints(species$range, nback, species$presence.points):
## generated random points = 0.9975 times requested number## Adding environmental data to species Iberolacerta monticola## Processing presence points...## Processing background points...## Warning in model.matrix.default(mt, mf, contrasts): non-list contrasts argument
## ignored## Fitting Models: 1 of 200
## Fitting Models: 2 of 200
## Fitting Models: 3 of 200
## Fitting Models: 4 of 200
## Fitting Models: 5 of 200
## Fitting Models: 6 of 200
## Fitting Models: 7 of 200
## Fitting Models: 8 of 200
## Fitting Models: 9 of 200
## Fitting Models: 10 of 200
## Fitting Models: 11 of 200
## Fitting Models: 12 of 200
## Fitting Models: 13 of 200
## Fitting Models: 14 of 200
## Fitting Models: 15 of 200
## Fitting Models: 16 of 200
## Fitting Models: 17 of 200
## Fitting Models: 18 of 200
## Fitting Models: 19 of 200
## Fitting Models: 20 of 200
## Fitting Models: 21 of 200
## Fitting Models: 22 of 200
## Fitting Models: 23 of 200
## Fitting Models: 24 of 200
## Fitting Models: 25 of 200
## Fitting Models: 26 of 200
## Fitting Models: 27 of 200
## Fitting Models: 28 of 200
## Fitting Models: 29 of 200
## Fitting Models: 30 of 200
## Fitting Models: 31 of 200
## Fitting Models: 32 of 200
## Fitting Models: 33 of 200
## Fitting Models: 34 of 200
## Fitting Models: 35 of 200
## Fitting Models: 36 of 200
## Fitting Models: 37 of 200
## Fitting Models: 38 of 200
## Fitting Models: 39 of 200
## Fitting Models: 40 of 200
## Fitting Models: 41 of 200
## Fitting Models: 42 of 200
## Fitting Models: 43 of 200
## Fitting Models: 44 of 200
## Fitting Models: 45 of 200
## Fitting Models: 46 of 200
## Fitting Models: 47 of 200
## Fitting Models: 48 of 200
## Fitting Models: 49 of 200
## Fitting Models: 50 of 200
## Fitting Models: 51 of 200
## Fitting Models: 52 of 200
## Fitting Models: 53 of 200
## Fitting Models: 54 of 200
## Fitting Models: 55 of 200
## Fitting Models: 56 of 200
## Fitting Models: 57 of 200
## Fitting Models: 58 of 200
## Fitting Models: 59 of 200
## Fitting Models: 60 of 200
## Fitting Models: 61 of 200
## Fitting Models: 62 of 200
## Fitting Models: 63 of 200
## Fitting Models: 64 of 200
## Fitting Models: 65 of 200
## Fitting Models: 66 of 200
## Fitting Models: 67 of 200
## Fitting Models: 68 of 200
## Fitting Models: 69 of 200
## Fitting Models: 70 of 200
## Fitting Models: 71 of 200
## Fitting Models: 72 of 200
## Fitting Models: 73 of 200
## Fitting Models: 74 of 200
## Fitting Models: 75 of 200
## Fitting Models: 76 of 200
## Fitting Models: 77 of 200
## Fitting Models: 78 of 200
## Fitting Models: 79 of 200
## Fitting Models: 80 of 200
## Fitting Models: 81 of 200
## Fitting Models: 82 of 200
## Fitting Models: 83 of 200
## Fitting Models: 84 of 200
## Fitting Models: 85 of 200
## Fitting Models: 86 of 200
## Fitting Models: 87 of 200
## Fitting Models: 88 of 200
## Fitting Models: 89 of 200
## Fitting Models: 90 of 200
## Fitting Models: 91 of 200
## Fitting Models: 92 of 200
## Fitting Models: 93 of 200
## Fitting Models: 94 of 200
## Fitting Models: 95 of 200
## Fitting Models: 96 of 200
## Fitting Models: 97 of 200
## Fitting Models: 98 of 200
## Fitting Models: 99 of 200
## Fitting Models: 100 of 200
## Fitting Models: 101 of 200
## Fitting Models: 102 of 200
## Fitting Models: 103 of 200
## Fitting Models: 104 of 200
## Fitting Models: 105 of 200
## Fitting Models: 106 of 200
## Fitting Models: 107 of 200
## Fitting Models: 108 of 200
## Fitting Models: 109 of 200
## Fitting Models: 110 of 200
## Fitting Models: 111 of 200
## Fitting Models: 112 of 200
## Fitting Models: 113 of 200
## Fitting Models: 114 of 200
## Fitting Models: 115 of 200
## Fitting Models: 116 of 200
## Fitting Models: 117 of 200
## Fitting Models: 118 of 200
## Fitting Models: 119 of 200
## Fitting Models: 120 of 200
## Fitting Models: 121 of 200
## Fitting Models: 122 of 200
## Fitting Models: 123 of 200
## Fitting Models: 124 of 200
## Fitting Models: 125 of 200
## Fitting Models: 126 of 200
## Fitting Models: 127 of 200
## Fitting Models: 128 of 200
## Fitting Models: 129 of 200
## Fitting Models: 130 of 200
## Fitting Models: 131 of 200
## Fitting Models: 132 of 200
## Fitting Models: 133 of 200
## Fitting Models: 134 of 200
## Fitting Models: 135 of 200
## Fitting Models: 136 of 200
## Fitting Models: 137 of 200
## Fitting Models: 138 of 200
## Fitting Models: 139 of 200
## Fitting Models: 140 of 200
## Fitting Models: 141 of 200
## Fitting Models: 142 of 200
## Fitting Models: 143 of 200
## Fitting Models: 144 of 200
## Fitting Models: 145 of 200
## Fitting Models: 146 of 200
## Fitting Models: 147 of 200
## Fitting Models: 148 of 200
## Fitting Models: 149 of 200
## Fitting Models: 150 of 200
## Fitting Models: 151 of 200
## Fitting Models: 152 of 200
## Fitting Models: 153 of 200
## Fitting Models: 154 of 200
## Fitting Models: 155 of 200
## Fitting Models: 156 of 200
## Fitting Models: 157 of 200
## Fitting Models: 158 of 200
## Fitting Models: 159 of 200
## Fitting Models: 160 of 200
## Fitting Models: 161 of 200
## Fitting Models: 162 of 200
## Fitting Models: 163 of 200
## Fitting Models: 164 of 200
## Fitting Models: 165 of 200
## Fitting Models: 166 of 200
## Fitting Models: 167 of 200
## Fitting Models: 168 of 200
## Fitting Models: 169 of 200
## Fitting Models: 170 of 200
## Fitting Models: 171 of 200
## Fitting Models: 172 of 200
## Fitting Models: 173 of 200
## Fitting Models: 174 of 200
## Fitting Models: 175 of 200
## Fitting Models: 176 of 200
## Fitting Models: 177 of 200
## Fitting Models: 178 of 200
## Fitting Models: 179 of 200
## Fitting Models: 180 of 200
## Fitting Models: 181 of 200
## Fitting Models: 182 of 200
## Fitting Models: 183 of 200
## Fitting Models: 184 of 200
## Fitting Models: 185 of 200
## Fitting Models: 186 of 200
## Fitting Models: 187 of 200
## Fitting Models: 188 of 200
## Fitting Models: 189 of 200
## Fitting Models: 190 of 200
## Fitting Models: 191 of 200
## Fitting Models: 192 of 200
## Fitting Models: 193 of 200
## Fitting Models: 194 of 200
## Fitting Models: 195 of 200
## Fitting Models: 196 of 200
## Fitting Models: 197 of 200
## Fitting Models: 198 of 200
## Fitting Models: 199 of 200
## Fitting Models: 200 of 200## Warning in model.matrix.default(mt, mf, contrasts): non-list contrasts argument
## ignored## Warning in model.matrix.default(mt, mf, contrasts): non-list contrasts argument
## ignored
## Warning in model.matrix.default(mt, mf, contrasts): non-list contrasts argument
## ignored
## Warning in model.matrix.default(mt, mf, contrasts): non-list contrasts argument
## ignored
## Warning in model.matrix.default(mt, mf, contrasts): non-list contrasts argument
## ignored
## Warning in model.matrix.default(mt, mf, contrasts): non-list contrasts argument
## ignored##
##
## Drawing background from species background points.## Warning in model.matrix.default(mt, mf, contrasts): non-list contrasts argument
## ignored
## Warning in model.matrix.default(mt, mf, contrasts): non-list contrasts argument
## ignored
## Warning in model.matrix.default(mt, mf, contrasts): non-list contrasts argument
## ignored
## Warning in model.matrix.default(mt, mf, contrasts): non-list contrasts argument
## ignored##
##
## Drawing background from species background points.## Warning in model.matrix.default(mt, mf, contrasts): non-list contrasts argument
## ignored
## Warning in model.matrix.default(mt, mf, contrasts): non-list contrasts argument
## ignored
## Warning in model.matrix.default(mt, mf, contrasts): non-list contrasts argument
## ignored
## Warning in model.matrix.default(mt, mf, contrasts): non-list contrasts argument
## ignored
## Warning in model.matrix.default(mt, mf, contrasts): non-list contrasts argument
## ignoredmonticola.ppm##
##
## Formula: ~bio3 + bio8 + bio11
## <environment: 0x7fcc668fd620>
##
##
## Data table (top ten lines):
##
## Longitude Latitude bio1 bio2 bio3 bio4 bio5 bio6 bio7 bio8 bio9 bio10 bio11 bio12 bio13 bio14 bio15 bio16 bio17 bio18 bio19 presence wt
## --- ---------- --------- ----- ----- ----- ----- ----- ----- ----- ----- ----- ------ ------ ------ ------ ------ ------ ------ ------ ------ ------ --------- ----------
## 2 -7.188726 43.43051 128 84 41 4080 239 38 201 83 180 182 78 959 124 33 34 345 126 142 310 1 0.0092593
## 3 -6.703289 43.46805 121 88 41 4237 239 28 211 97 175 178 70 920 121 36 31 320 135 148 279 1 0.0039683
## 4 -6.504768 43.51037 134 85 42 3987 243 45 198 113 185 187 86 860 116 35 30 299 130 145 248 1 0.0092593
## 5 -6.556106 43.52801 134 85 42 3987 243 45 198 113 185 187 86 860 116 35 30 299 130 145 248 1 0.0092593
## 6 -6.728039 43.55888 137 83 42 3957 244 48 196 115 187 190 89 852 114 33 31 298 124 139 253 1 0.0055556
## 7 -6.823801 43.39811 121 88 41 4237 239 28 211 97 175 178 70 920 121 36 31 320 135 148 279 1 0.0039683
## 8 -6.732560 43.54083 137 83 42 3957 244 48 196 115 187 190 89 852 114 33 31 298 124 139 253 1 0.0055556
## 9 -6.650406 43.40437 117 90 42 4292 237 23 214 93 172 174 65 926 123 38 30 320 140 153 274 1 0.0069444
## 10 -6.541960 43.45676 117 90 42 4292 237 23 214 93 172 174 65 926 123 38 30 320 140 153 274 1 0.0069444
## 11 -6.910225 43.29794 117 91 41 4472 240 20 220 67 175 176 63 942 122 34 33 334 129 142 300 1 0.0034722
##
##
## Model: Length Class Mode
## betas 804 -none- numeric
## lambdas 201 -none- numeric
## likelihoods 201 -none- numeric
## pen.likelihoods 201 -none- numeric
## lambda 1 -none- numeric
## beta 4 -none- numeric
## mu 645 -none- numeric
## likelihood 1 -none- numeric
## criterion 1 -none- character
## family 1 -none- character
## gamma 1 -none- numeric
## alpha 1 -none- numeric
## init.coef 1 -none- logical
## criterion.matrix 5 data.frame list
## data 2580 -none- numeric
## pt.interactions 1 -none- logical
## wt 645 -none- numeric
## pres 645 -none- numeric
## x 0 -none- NULL
## y 0 -none- NULL
## r 1 -none- logical
## call 5 -none- call
## formula 2 formula call
## s.means 3 -none- numeric
## s.sds 3 -none- numeric
## cv.group 645 -none- numeric
## n.blocks 1 -none- logical
##
##
## Environment space model fit (training data): class : ModelEvaluation
## n presences : 246
## n absences : 10000
## AUC : 0.6614008
## cor : -0.02189001
## max TPR+TNR at : 17.99064
##
##
## Proportion of data wittheld for model fitting: 0.3
##
## Model fit (test data): class : ModelEvaluation
## n presences : 106
## n absences : 399
## AUC : 0.8606303
## cor : 0.5689826
## max TPR+TNR at : 22.86892
##
##
## Environment space model fit (test data): class : ModelEvaluation
## n presences : 106
## n absences : 10000
## AUC : 0.7373255
## cor : 0.004216513
## max TPR+TNR at : 17.99064
##
##
## Suitability:
## class : RasterLayer
## dimensions : 60, 84, 5040 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -10, 4, 35, 45 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : layer
## values : 0.0002005848, 28.77812 (min, max)
##
##
##
## Notes:
monticola.ppm$response.plots## $bio3
##
## $bio8
##
## $bio11
visualize.enm(monticola.ppm, spain.worldclim, layers = c("bio3", "bio8"))## Warning in model.matrix.default(mt, mf, contrasts): non-list contrasts argument
## ignored## $background.plot## Warning: Removed 396 rows containing missing values (geom_raster).
##
## $suit.plot
One final option is to use Blonder et al.’s hypervolume package. At present these are implemented as enmtools.hypervolume objects, rather than as enmtools.model objects. They don’t have as many cool options or visualizations, as many of the options for models don’t make much sense here. It’s worth knowing that you can build them, though, as we will eventually be using them for some of the hypothesis tests.
enmtools.hypervolume(monticola, spain.worldclim[[c(1,8,12,17)]])##
## Building tree...
## done.
## Ball query...
##
## done.
## Requested probability quantile 0.950000, obtained 0.947436 - setting threshold value 0.000025.
## For a closer match, you can increase num.thresholds in hypervolume_threshold.
## Retaining 172/1726 hypervolume random points for comparison with 5040 test points.## $hv
## ***** Object of class Hypervolume *****
## Name: Iberolacerta monticola
## Method: Gaussian kernel density estimate
## Number of data points (after weighting): 352
## Dimensionality: 4
## Volume: 14.850743
## Random point density: 116.223139
## Number of random points: 1726
## Random point values:
## min: 0.000
## mean: 0.000
## median: 0.000
## max:0.002
## Parameters:
## kde.bandwidth: 0.2444054 0.2731535 0.2140918 0.1427126
## samples.per.point: 10
## sd.count: 3
## quantile.requested: 0.95
## quantile.requested.type: probability
##
## $suitability
## class : RasterLayer
## dimensions : 60, 84, 5040 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -10, 4, 35, 45 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : bio1
## values : 0, 0.0002983051 (min, max)
##
##
## $species.name
## [1] "Iberolacerta monticola"
##
## $analysis.df
## Longitude Latitude
## 1 -6.699292 43.45172
## 2 -7.188726 43.43051
## 3 -6.703289 43.46805
## 4 -6.504768 43.51037
## 5 -6.556106 43.52801
## 6 -6.728039 43.55888
## 7 -6.823801 43.39811
## 8 -6.732560 43.54083
## 9 -6.650406 43.40437
## 10 -6.541960 43.45676
## 11 -6.910225 43.29794
## 12 -6.659564 43.48790
## 13 -6.817537 43.28005
## 14 -6.831939 43.38955
## 15 -6.793121 43.45933
## 16 -6.499937 43.47067
## 17 -6.619466 43.42050
## 18 -6.724369 43.52515
## 19 -7.657592 40.33878
## 20 -6.660907 43.55608
## 21 -6.933462 43.48396
## 22 -6.687861 43.43833
## 23 -6.951364 43.59288
## 24 -6.688261 43.39621
## 25 -6.704228 43.36865
## 26 -6.619389 43.51894
## 27 -8.150035 43.34445
## 28 -6.618920 43.55074
## 29 -6.770702 43.50165
## 30 -7.754594 40.34196
## 31 -6.770915 43.49353
## 32 -4.896618 43.24043
## 33 -7.654091 40.39647
## 34 -7.922660 43.77674
## 35 -6.746105 43.48776
## 36 -5.975720 43.37362
## 37 -6.958813 43.35228
## 38 -6.604129 43.21463
## 39 -6.674231 43.23169
## 40 -7.049430 42.79107
## 41 -7.154147 42.62575
## 42 -7.128125 42.77214
## 43 -6.231050 43.19318
## 44 -7.056323 43.39149
## 45 -5.621583 43.01901
## 46 -7.982539 43.65810
## 47 -7.810377 43.63642
## 48 -7.824902 43.61708
## 49 -6.904581 43.29112
## 50 -6.827700 42.96513
## 51 -8.040703 43.54231
## 52 -8.008715 43.51093
## 53 -6.931485 43.26939
## 54 -6.763920 43.59477
## 55 -6.927234 43.24961
## 56 -6.180002 43.19246
## 57 -7.161566 42.67164
## 58 -7.740096 43.75711
## 59 -7.480695 40.32510
## 60 -7.729511 40.39481
## 61 -6.264794 43.04339
## 62 -6.276012 43.11746
## 63 -6.964268 43.25487
## 64 -6.745783 43.28473
## 65 -6.898205 43.53368
## 66 -6.139212 43.05136
## 67 -6.142102 42.93155
## 68 -5.171215 43.06957
## 69 -6.036635 43.02531
## 70 -6.651737 43.47214
## 71 -6.690762 43.45555
## 72 -7.473340 43.78935
## 73 -7.411094 43.76633
## 74 -6.575039 42.91070
## 75 -8.180498 43.38845
## 76 -8.063401 43.24602
## 77 -8.053342 43.23926
## 78 -6.834718 42.89244
## 79 -4.878881 43.29832
## 80 -6.053402 43.15157
## 81 -6.050337 43.07146
## 82 -6.132320 43.00186
## 83 -8.144452 43.39255
## 84 -7.787817 40.24217
## 85 -4.686575 43.01294
## 86 -7.150816 43.17367
## 87 -5.132756 43.49572
## 88 -7.787378 40.39362
## 89 -7.340277 42.74906
## 90 -4.812590 43.36524
## 91 -8.110797 43.70124
## 92 -4.941888 43.35310
## 93 -6.971264 43.50353
## 94 -7.645674 40.36543
## 95 -7.621731 40.34170
## 96 -7.642539 40.36317
## 97 -6.803279 43.30922
## 98 -6.968332 43.34016
## 99 -6.961694 43.34072
## 100 -4.780000 43.02000
## 101 -4.780000 43.11000
## 102 -4.780000 43.20000
## 103 -4.790000 43.29000
## 104 -4.910000 43.11000
## 105 -4.910000 43.20000
## 106 -5.030000 43.11000
## 107 -5.030000 43.20000
## 108 -5.150000 43.11000
## 109 -5.270000 43.11000
## 110 -5.390000 43.01000
## 111 -5.400000 43.10000
## 112 -5.520000 43.01000
## 113 -5.520000 43.10000
## 114 -5.760000 43.00000
## 115 -5.880000 43.00000
## 116 -6.110000 43.09000
## 117 -6.120000 43.00000
## 118 -6.240000 43.00000
## 119 -6.240000 43.09000
## 120 -6.360000 43.10000
## 121 -6.520000 42.20000
## 122 -6.610000 42.92000
## 123 -6.640000 42.20000
## 124 -6.760000 42.21000
## 125 -6.760000 42.30000
## 126 -6.850000 43.02000
## 127 -6.850000 43.20000
## 128 -6.860000 42.75000
## 129 -6.860000 42.84000
## 130 -6.860000 42.93000
## 131 -6.870000 42.39000
## 132 -6.880000 42.12000
## 133 -6.970000 43.20000
## 134 -6.980000 42.75000
## 135 -6.980000 43.02000
## 136 -6.990000 42.57000
## 137 -6.990000 42.66000
## 138 -7.010000 41.94000
## 139 -7.080000 43.47000
## 140 -7.090000 43.20000
## 141 -7.090000 43.29000
## 142 -7.090000 43.38000
## 143 -7.110000 42.57000
## 144 -7.210000 43.38000
## 145 -7.240000 42.30000
## 146 -7.700000 43.48000
## 147 -7.700000 43.66000
## 148 -7.700000 43.75000
## 149 -7.830000 43.48000
## 150 -7.950000 43.30000
## 151 -7.950000 43.39000
## 152 -7.960000 43.03000
## 153 -7.820000 43.57000
## 154 -7.820000 43.66000
## 155 -7.820000 43.75000
## 156 -7.940000 43.75000
## 157 -7.950000 43.48000
## 158 -7.950000 43.57000
## 159 -7.950000 43.66000
## 160 -8.070000 43.39000
## 161 -8.070000 43.48000
## 162 -8.070000 43.66000
## 163 -8.070000 43.75000
## 164 -8.080000 43.03000
## 165 -8.080000 43.21000
## 166 -8.080000 43.30000
## 167 -8.200000 43.31000
## 168 -8.200000 43.40000
## 169 -8.320000 43.31000
## 170 -9.060000 42.86000
## 171 -5.150000 43.02000
## 172 -5.270000 43.02000
## 173 -5.390000 42.92000
## 174 -5.510000 42.83000
## 175 -5.630000 42.92000
## 176 -6.390000 42.29000
## 177 -6.610000 42.83000
## 178 -6.630000 42.29000
## 179 -6.740000 42.75000
## 180 -6.740000 42.84000
## 181 -6.980000 42.84000
## 182 -7.090000 43.11000
## 183 -7.100000 42.84000
## 184 -7.100000 42.93000
## 185 -7.100000 43.02000
## 186 -7.110000 42.66000
## 187 -7.210000 43.29000
## 188 -7.330000 43.39000
## 189 -7.330000 43.48000
## 190 -7.350000 42.67000
## 191 -7.450000 43.48000
## 192 -7.570000 43.66000
## 193 -7.580000 43.48000
## 194 -7.240000 42.21000
## 195 -7.360000 42.22000
## 196 -7.490000 42.22000
## 197 -4.910000 43.29000
## 198 -5.030000 43.29000
## 199 -5.160000 43.47000
## 200 -5.280000 43.20000
## 201 -5.290000 43.47000
## 202 -5.400000 43.19000
## 203 -5.520000 43.19000
## 204 -5.890000 43.09000
## 205 -5.890000 43.18000
## 206 -6.110000 43.18000
## 207 -6.480000 43.10000
## 208 -6.600000 43.28000
## 209 -6.610000 43.01000
## 210 -6.710000 43.47000
## 211 -6.720000 43.20000
## 212 -6.720000 43.29000
## 213 -6.730000 43.02000
## 214 -6.730000 43.11000
## 215 -6.840000 43.38000
## 216 -6.840000 43.47000
## 217 -6.960000 43.38000
## 218 -6.960000 43.47000
## 219 -6.960000 43.56000
## 220 -6.970000 43.29000
## 221 -6.760000 42.12000
## 222 -5.820000 43.14000
## 223 -7.613616 40.32185
## 224 -4.910000 43.15000
## 225 -7.363881 43.42967
## 226 -6.430000 42.88000
## 227 -6.190000 42.85000
## 228 -6.550000 43.18000
## 229 -6.660000 42.11000
## 230 -7.616667 40.31667
## 231 -7.638880 40.36111
## 232 -6.440000 42.26000
## 233 -6.730000 42.82000
## 234 -6.620000 42.77000
## 235 -6.580000 42.27000
## 236 -5.530000 42.89000
## 237 -0.020000 42.59000
## 238 -2.450000 42.46000
## 239 -7.630000 40.33000
## 240 -8.080000 43.03000
## 241 -8.320000 43.31000
## 242 -9.060000 42.86000
## 243 -8.200000 43.31000
## 244 -8.080000 43.21000
## 245 -7.960000 43.03000
## 246 -7.950000 43.30000
## 247 -8.080000 43.30000
## 248 -7.090000 43.29000
## 249 -6.970000 43.20000
## 250 -7.950000 43.39000
## 251 -8.070000 43.39000
## 252 -8.200000 43.40000
## 253 -7.090000 43.20000
## 254 -7.090000 43.11000
## 255 -6.970000 43.29000
## 256 -7.240000 42.21000
## 257 -6.630000 42.29000
## 258 -7.210000 43.29000
## 259 -7.490000 42.22000
## 260 -7.350000 42.67000
## 261 -6.870000 42.39000
## 262 -7.240000 42.30000
## 263 -7.360000 42.22000
## 264 -6.760000 42.30000
## 265 -6.760000 42.21000
## 266 -6.760000 42.12000
## 267 -7.010000 41.94000
## 268 -6.640000 42.20000
## 269 -6.880000 42.12000
## 270 -5.630000 42.92000
## 271 -4.830000 40.50000
## 272 -5.510000 42.83000
## 273 -5.510000 42.92000
## 274 -4.710000 40.32000
## 275 -5.650000 40.31000
## 276 -6.340000 43.46000
## 277 -5.520000 43.01000
## 278 -5.520000 43.19000
## 279 -5.970000 42.99000
## 280 -5.770000 40.30000
## 281 -5.520000 43.10000
## 282 -5.180000 40.41000
## 283 -5.970000 43.09000
## 284 -4.820000 40.32000
## 285 -6.610000 42.92000
## 286 -7.200000 43.56000
## 287 -6.480000 43.10000
## 288 -6.720000 43.20000
## 289 -6.360000 43.01000
## 290 -6.720000 43.29000
## 291 -6.610000 42.83000
## 292 -6.610000 43.01000
## 293 -6.370000 42.92000
## 294 -6.710000 43.47000
## 295 -7.580000 43.57000
## 296 -6.840000 43.38000
## 297 -6.850000 43.20000
## 298 -6.960000 43.56000
## 299 -7.330000 43.39000
## 300 -7.210000 43.38000
## 301 -7.080000 43.47000
## 302 -6.730000 43.11000
## 303 -6.960000 43.47000
## 304 -6.520000 42.20000
## 305 -6.840000 43.47000
## 306 -6.110000 40.57000
## 307 -7.090000 43.38000
## 308 -6.390000 42.29000
## 309 -5.990000 43.27000
## 310 -6.360000 43.10000
## 311 -6.600000 43.28000
## 312 -6.350000 43.19000
## 313 -6.960000 43.38000
## 314 -6.010000 43.13000
## 315 -6.020000 42.99000
## 316 -4.910000 43.11000
## 317 -5.150000 43.11000
## 318 -4.910000 43.20000
## 319 0.130000 42.64000
## 320 -5.390000 42.92000
## 321 -5.280000 43.20000
## 322 -3.890000 40.87000
## 323 -5.030000 43.11000
## 324 -5.030000 43.20000
## 325 -4.780000 43.11000
## 326 -5.400000 43.10000
## 327 -4.910000 43.29000
## 328 -5.390000 43.01000
## 329 -4.130000 40.78000
## 330 -3.890000 40.96000
## 331 -3.770000 40.78000
## 332 -5.030000 43.29000
## 333 -4.790000 43.29000
## 334 -5.270000 43.11000
## 335 -4.660000 43.21000
## 336 -3.890000 40.69000
## 337 -4.590000 40.32000
## 338 -4.780000 43.20000
## 339 -5.160000 43.47000
## 340 -3.770000 41.05000
## 341 -5.290000 43.47000
## 342 -5.400000 43.19000
## 343 -4.010000 40.69000
## 344 -4.010000 40.78000
## 345 -3.770000 40.96000
## 346 -4.650000 42.94000
## 347 -5.280000 43.38000
## 348 -3.890000 40.78000
## 349 -4.120000 40.60000
## 350 -5.160000 43.38000
## 351 -3.650000 41.05000
## 352 -4.010000 40.87000
##
## attr(,"class")
## [1] "enmtools.hypervolume"Projections across space and time
Often we have applications where we want to take a model made in one time and place and project it to the distribution of environments in another time and place. This is relatively easy to do in ENMTools! For instance we can extrapolate our model across the entire world.
mx.pred <- predict(monticola.mx, all.worldclim)
mx.pred$suitability## Warning: Raster pixels are placed at uneven vertical intervals and will be
## shifted. Consider using geom_tile() instead.
It’s a bit crazy to think that these lizards are going to live at the north pole, so let’s look at the built-in plots that visualize the environments where the model was trained vs. where it’s being projected.
mx.pred$threespace.plot
We can look at the suitability projected in the future as well. Note: this involves downloading new layers so this bit of code may take a while to run.
spain.future <- raster::getData('CMIP5', var = "bio", res = 10, rcp = 85, model = 'AC', year = 70)
spain.extent <- extent(-10, 4, 35, 45)
spain.future <- crop(x = spain.future, y = spain.extent)
names(spain.future) <- names(spain.worldclim)
future.pred <- predict(monticola.mx, spain.future)
future.pred$suitability
future.pred$threespace.plot
Note that you can call predict on the model objects stored inside ENMTools models as well, but you often need to change the order. For instance to predict using a maxent model, you pass in the environmental rasters first and the model second. This is just a dismo thing to be aware of, there’s no particular logic to either ordering that you need to understand.
plot(predict(all.worldclim, monticola.mx$model))